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Abstract. The capture and gradual inspiral of stellar mass objects by a massive 

black hole at the centre of a galaxy has been proposed as one of the most promising 

source of gravitational radiation to be detected by LISA. Unfortunately rate estimates 

for this process suffer from many uncertainties. Here we report on the use of our newly 

^. . developed Monte Carlo stellar dynamics code to tackle this problem. We present 

C*"") | results from simple galactic nuclei models that demonstrate the high potential of our 

O^ ■ approach and point out the aspects of the problem where an improved treatment seems 

■^^ , desirable. 
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PACS numbers: 98.62.Js, 02.70.Lq, 95.85.Sz, 98.10. +z, 97.60.-s 

1. Introduction 

In the past decade, the harvest of observational evidences for the presence of massive 
dark objects in the centre of most bright galaxies has become very impressive. It is 
generally accepted that these mass concentrations are most likely massive black holes 
(MBHs) with detected masses ranging between 10 6 and 1O 1O M . For all but two 
galaxies, however, today's observations in the electro-magnetic spectrum lack resolution 
to rule out concurrent interpretations concerning the nature of the dark object |], [|. 

The motion of a "test particle" (with mass -C Mmbh) on a relativistic orbit, 
monitored by the emitted gravitational waves (GW), is an ideal way to get much more 
precise information from the immediate vicinity of the putative black hole. By detecting 
such signals and comparing them with general relativity predictions, we could, at the 
same time, test Einstein's theory in the strong field regime, confirm the existence of 
MBHs (in a statistical sense, if not in particular galaxies) and get precise determinations 
of their mass and spin || . 

However, two main classes of theoretical difficulties have to be worked around. First, 
the computation of the trajectory of a particle around a Kerr BH without resorting to 
strong simplifying assumptions, is a formidable, still unsolved, problem. The second 
difficulty is the prediction of the number of potential sources, i.e. the rate with which 
stars get so deeply bound to a central MBH that their subsequent orbital evolution is 
driven by GW-emission. 
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The range of MBH masses that can be probed by LISA is 3 x 10 5 - 3 x 10 7 M Q §]. 
Only for the highest Mmbh values can MS stars be captured on relativistic orbits without 
being teared off by tidal forces so we expect most capture events to feature a compact 
remnant. It can be brought onto relativistic orbits through 2-body relaxation or non 
disruptive collisions with other stars if the velocity dispersion exceeds V*, an unlikely 
situation for compact stars. Repeated impacts on an accretion disk could also channel 
stars onto relativistic orbits @j, but this process is probably ineffective for compact 
remnants as they would not sweep much disk's material at each crossing. 

The rate of captures on orbits with strong emission of gravitational radiation 
(hereafter "GR-captures" ) has been estimated in three previous works f| || [/J. 
Unfortunately, these authors use different galactic nuclei models, assumptions and 
computation methods. . . and quite expectedly get discordant capture rates, ranging 
from 2 x 10~ 8 to 10 yr -1 . Besides the determination of such a rate, important issues 
raised by these works are the following. To what extent is the mass-segregation process 
efficient at concentrating stellar remnants towards the central regions? Is the GR- 
capture rate sustainable for many Gyrs? How does the stellar nucleus evolve on the 
long run in response to this and other stellar dynamical processes (relaxation, tidal 
disruptions, collisions, stellar evolution)? 

2. The Monte Carlo cluster evolution code 

2.1. Overview of the code 

We recently developed a new computer code to follow the evolution of a central star 
cluster surrounding a MBH over time scales as long as 10 10 years. This tool, which 
we describe in ||, is the backbone of our work aimed at a better modelization and 
understanding of the influence of the MBH on the stellar system. 

The numerical scheme we chose is a Monte Carlo (MC) method based on the 
pioneering work of Henon M. This class of programs proved very successful in the 



realm of globular clusters studies [10, 11, 12]. In its basic form, the method is more 



or less equivalent to resolving the Fokker-Planck (FP) equation, but with no implicit 
integration of the diffusion coefficients. The relaxation-driven secular evolution of a 
stellar system is simulated without wasting CPU time at integrating orbits. Compared 
to direct FP resolutions, the main advantages of the MC method are (1) its simplicity, 

(2) the fact that any stellar mass function and any velocity distribution (both variable 
in time and position) are naturally coped with, as is the cluster's self-gravity, and 

(3) the ease with which extra physical processes can be included, even those that are 
not continuous, like collisions. Although they are also particle-based, MC simulations 
are enormously faster than relaxational A-body integrations whose computing time 
increases like Np~l while the scaling is almost linear for our MC code. Furthermore, as 
each particle represents a spherical shell of stars with same orbital and stellar properties, 
the number of simulated stars can be made arbitrarily high. Of course, MC methods 
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also suffer from a few drawbacks. On the one hand, they produce much noisier results 
than FP methods. On the other hand, they rely on more simplifying assumptions than 
iV-body schemes; the strongest of them are spherical symmetry, constant dynamical 
equilibrium and treatment of relaxation as the integrated effect of distant, uncorrelated, 
2-body encounters. 

As far as we know, our code is the first implementation of the Henon's scheme 
devoted to MBH-hosting dense galactic nuclei. We spent much time to include stellar 
collisions (between MS stars, for the time being) as realistically as possible |13, |14 |. From 



more than 12 000 MS-MS collision simulations carried out with a Smoothed Particle 
Hydrodynamics code (SPH, [[0|), we built a grid of collisional outcomes from which the 
result of any collision is interpolated when it occurs during a cluster evolution simulation. 
We also included "loss-cone" effects (TBI to account for the tidal disruptions of stars in 



the Roche zone around the MBH. Finally, for this LISA symposium, we adapted this 
loss-cone routine in order to detect GR-capture events. 

2.2. Computing capture rates with the Monte Carlo code 

Following Sigurdsson & Rees ||, we consider that a star is "captured" by the MBH 
if it gets on an orbit whose shrinkage time by GW-emission, Tqw, is shorter than the 
time 2-body relaxation would take to modify it substantially, T mo d. The most likely way 
for a star to get sufficiently close to the MBH to be captured is to experience 2-body 
encounters that make its orbit nearly radial. For such a highly eccentric orbit, T mot j is 
the time taken by relaxation to introduce a relative change in the pericentre distance of 
order 1. Then, T mo d — (1 — e)T rc i ax <C T rc i ax where e is the eccentricity. Typical values 
of 1 — e for captured particles are 10~ 7 — 10~ 3 . 

Using time steps 5t as small as a fraction of T mod would completely jeopardize the 
efficiency of the MC scheme for which 5t ~ 0.001 — 0.01T re i ax is otherwise sufficiently 
small. For a star with fixed \V\ at a given distance from the centre, the capture orbits 
form a thin conical bundle with aperture angle 6*grc <C it. If we impose too large a 5t, V 
jumps over this tiny "loss cone" and most captures are missed. This problem is closely 
reminiscent to the detection of tidal disruptions. Thus we solved it the same way. We 
keep time steps that are a (small) fraction of T re i ax (and/or T co u) but we "over-sample" 
each of them by simulating the random walk of the tip of V on a sphere during St. 
We simply piggy-backed the procedure for GR-captures detections on the routine for 
tidal disruptions. Hence, as a star would be tidally destroyed only if it stays on a loss 
cone orbit until next pericentre passage, the time resolution of the random walk is the 
orbital period T orb . This is not well adapted to GR-captures for which T orb has not 
special meaning. However, entries in the GR-capture loss cone are correctly simulated 
as a diffusion process as long as individual steps in the random walk are smaller than 
#gro But this is always true because, by definition, the R.M.S. diffusion angle during 
Tgw is #grc and Tqw is much longer than T or t,. 

We think that the main intrinsic limitation in our approach is to treat relaxation 
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with small angle approximation. Indeed, close gravitational encounters with scattering 
angles of order n, although they only contribute a fraction ln(A) -1 < 0.1 to the overall 
relaxation, can dominate the rate of captures ||. Another limitation is of statistical 
nature. To ensure energy conservation, each particle stands for a fixed number of stars. 
Consequently there are only 1300 SBH particles out of 2 x 10 6 . 

3. Simulations of simple models of galactic nuclei 

3.1. Initial model & included physics 

To explore the potential of our MC code in predicting GR-capture rates, we run a few 
simulations that are variations around our "standard" galactic nucleus model. The 
stellar cluster is set as a Wq = 8 King model, with a core radius of 0.47 pc and 3.6 x 10 8 
stars. As stellar evolution is not simulated, we have to include an evolved stellar 
population, containing compact remnants, from the beginning of the computations. It 
is prepared according to the prescription of Miralda-Escude & Gould |7j . The number 
fractions of white dwarfs (WDs), neutron stars (NSs) and SBHs are 0.056, 3.5 x 10~ 3 
and 6.5 x 10~ 4 , and their individual masses 0.6, 1.4 and 7.0 M , respectively. We do not 
include giant stars. An initial "seed" black hole (MBn/M c \ nst ~ 10~ 4 ) is set at the centre. 
It accretes all gas released by stellar collisions and tidal disruptions with no delay. GR- 
captured stars are also assumed to be immediately swallowed by the BH. Direct plunges 
through the horizon are detected. Collisions between MS stars are treated realistically 
thank to our SPH grid but tidal disruptions are assumed to be always total. 

This model is not meant to be highly realistic. It is probably too dense and massive 
to represent a real galactic nucleus. To bracket the density value at the Galactic centre 
|j~7|l , we also used another model with a number of stars reduced by a factor 10. A 



better representation for a galactic nucleus would require a power law density profile 
at large radii (p oc R~ a with a < 2 typically) instead of a steep cut-off. This point 
is problematical as it imposes to put large amounts of particles at large distances and 
consequently reduce the resolution near the MBH where all the action takes place. 

3.2. Simulation results 

Some aspects of the evolution of the "standard" model are shown in figure [3] while 
figure H illustrates the evolution of the lighter nucleus. Both simulations were realized 
with 2 x 10 6 particles. 

We see that relaxational segregation of SBHs towards the centre occurs quickly. An 
early broad peak in the capture rate for this species ensues. Asbh then steadily decreases 
due to exhaustion of SBHs in the central regions. At late times, the capture rate is 
dominated by WDs, NSs, and, for dense clusters, low-mass MS stars ((M*) ~ 0.2 M ). 
However, we noted that for many of these MS stars, we have Tqw ^ 10 9 years even 
though the capture condition (Tqw < T moA ) is obeyed. Further investigations about 
these cases are called for. Typical remnant capture rates at T = 1 — 2 x 10 10 years are 
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Figure 1. Evolution of the "standard" model (see text). Panel (a) shows the rate 
of captures by emission of gravitational radiation for various stellar types. Panel (b) 
depicts the evolution of a few Lagrangian radii. Panel (c) shows the contribution of 
various processes to the growth of the central BH. 
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Figure 2. Same as figure p] for a nucleus model which is 10 times less massive. 

~ 10 _5 yr _1 and ~5x 10 _7 yr _1 for the standard and "light" models. A nucleus with 
10 times the mass of the standard model yields ~ 10 _4 yr _1 . 

A massive BH (M > 10 6 M Q ) can be grown from a seed if the initial stellar density 
is high enough. However, the stellar density after a Hubble time remains uncomfortably 
high (~ 10 7 M Q pc~ 3 at 0.1 pc). Tidal disruptions are the major contributors to the BH's 
mass with stellar collisions playing only a minor role. In the lighter model, GR-captures 
of SBHs dominates the early growth of the central BH with tidal disruptions catching 
up later on. This is probably connected to the stronger relaxation in this model. 



3.3. Future improvements 

There are many ways our simulations could be improved in the future. Let's mention a 
few of them that are relevant for the problem of GR-captures. 

• The lack of resolution for rare species, such as SBHs, and the problem of extending 
the initial model to larger radii could be solved if each particle could represent 
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a different number of stars. We presently do not know how to achieve this. 
Alternatively, a massive increase in the number of particles would solve these 
difficulties. A simulation with 10 7 particles would require ~ 1 month of CPU 
time, but we run short of computer memory. To go beyond this number, a parallel 
version of the code should be developed. 

• Stellar evolution should be included. In taht case, we have to determine the fate 
of the gas lost by stars, and particularly the fraction accreted onto the central BH. 
For our IMF, more than 40 % of the ZAMS mass is lost from the stars in the first 
10 10 years, which amounts to 5.9 x 10 7 M Q for the standard model! Preliminary 
calculations including stellar evolution indicate that the GR-capture rate at late 
stages is not substantially altered when all the emitted gas is lost from the nucleus. 

• A better treatment of GR-capture events should be attempted. It could go beyond 
the simple Tqw < ^mod criterion and the assumption of immediate swallowing. 

• The contribution of large angle scatterings to GR-captures must be clarified and 
included in our simulations if important. 
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